function interpolation_comparison()
    % 原始数据点
    x = [0 1 4 9 16 25 36 49 64];
    y = [0 1 2 3  4  5  6  7  8];
    
    % 创建更密的点用于绘图
    x_plot = linspace(0, 64, 200);
    
    % 1. Newton插值
    y_newton = newton_interp(x, y, x_plot);
    
    % 2. 分段线性插值
    y_linear = interp1(x, y, x_plot, 'linear');
    
    % 3. 3次最小二乘拟合
    p = polyfit(x, y, 3);
    y_least_squares = polyval(p, x_plot);
    
    % 绘图
    figure('Position', [100, 100, 800, 600]);
    plot(x, y, 'ko', 'MarkerSize', 8, 'LineWidth', 1.5); % 原始数据点
    hold on;
    plot(x_plot, y_newton, 'r-', 'LineWidth', 1.5, 'DisplayName', 'Newton插值');
    plot(x_plot, y_linear, 'b--', 'LineWidth', 1.5, 'DisplayName', '分段线性插值');
    plot(x_plot, y_least_squares, 'g-.', 'LineWidth', 1.5, 'DisplayName', '3次最小二乘拟合');
    plot(x_plot, sqrt(x_plot), 'k:', 'LineWidth', 1.5, 'DisplayName', '原函数');
    
    grid on;
    legend('数据点', 'Newton插值', '分段线性插值', '3次最小二乘拟合', '原函数');
    title('f(x)=\surd{x} 的不同逼近方法比较');
    xlabel('x');
    ylabel('y');
    
    % 计算并输出误差
    fprintf('\n各方法在插值点处的最大误差：\n');
    err_newton = max(abs(newton_interp(x, y, x) - y));
    err_linear = max(abs(interp1(x, y, x, 'linear') - y));
    err_least = max(abs(polyval(p, x) - y));
    
    fprintf('Newton插值最大误差: %.6e\n', err_newton);
    fprintf('分段线性插值最大误差: %.6e\n', err_linear);
    fprintf('3次最小二乘拟合最大误差: %.6e\n', err_least);
end

function y_interp = newton_interp(x, y, x_interp)
    n = length(x);
    % 计算差商
    d = zeros(n, n);
    d(:,1) = y(:);
    
    for j = 2:n
        for i = 1:(n-j+1)
            d(i,j) = (d(i+1,j-1) - d(i,j-1))/(x(i+j-1) - x(i));
        end
    end
    
    % 计算插值结果
    y_interp = d(1,1) * ones(size(x_interp));
    for j = 2:n
        term = d(1,j);
        for k = 1:(j-1)
            term = term .* (x_interp - x(k));
        end
        y_interp = y_interp + term;
    end
end